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1 Introduction 

In recent years we have witnessed rapid progress in the understanding of the mathe- 
matical structure of the static correlation functions of Yang-Baxter integrable quantum 
systems. Most of this progress was obtained with the example of the XXZ spin-^ 
chain. For the XXZ chain a hidden Grassmann structure was identified in [5, 6] which 
made it possible to prove the complete factorization of the correlation functions under 
very general conditions [19], including the case of finite temperature and magnetic 
field [2,19]. 

At the outset of this new development were multiple integral representations for the 
density matrix of a finite chain segment [14, 17, 18,24] and the observation in [7] that 
these integrals factorize into sums over products of single integrals. With [3,9] it became 
apparent that the factorization is not a property of the ground state in the thermodynamic 
limit, but can be done for finite temperature and for the ground states of finite chains 
as well. This was part of the motivation for the research leading to [2, 5, 6, 19]. The 
multiple integral representations also served as the starting point for a direct calculation 
of the asymptotics of the ground state correlation functions of the XXZ chain in [23]. 

At the present stage of research it is an interesting question to which extend the 
results for the XXZ spin-^ chain can be generalized to other integrable models. The 
models closest to the spin- 1 XXZ chain are those with the same /^-matrix, notably 
the Bose gas and the Sine-Gordon model. For both of these, partial results could be 
obtained [20, 22, 28, 29] in certain scaling limits. Another class of models, which 
is closely related to the spin-j XXZ chain as well, is the class of its higher-spin 
generalizations constructed by means of the fusion procedure [26,27]. 

For the fused spin chains N. Kitanine constructed a multiple integral representation 
[21] for the ground state correlation functions. He observed that much of the necessary 
algebraic and combinatorial work can be carried over rather directly from the spin- ^ 
case [24]. But due to the different structure of the ground state, which is build up 
of strings of Bethe roots for the higher spin integrable chains, the rewriting of the 
combinatorial sums as integrals in the thermodynamic limit required some modification 
as compared to the spin-j case. As a result the number of integrals in Kitanine's formula 
is 2ms for the m-site density matrix of the spin-.? chain, and a subtle regularization 
determines the relative location of the integration contours. Unlike in the spin-j 
case his multiple integral formula for higher spins bears no obvious similarity with the 
formulae obtained within the ^-vertex operator approach [8, 16]. For simplicity Kitanine 
concentrated on the isotropic (or XXX-) case, and he did not include a magnetic field. 
The generalization of his work to the XXZ-case (without magnetic field) was recently 
obtained in [10]. 

It is the aim of this work to extend Kitanine's result, exemplarily in the simplest 
case of the isotropic spin-1 chain, to finite temperatures. We shall also include a 
magnetic field into the calculation. Again fusion allows us to start with spin- \ and 
to use the algebraic and combinatorial results of [12,24]. Then, as we shall see, the 
crucial problem is the analytic part of the calculation, where the combinatorial sums 
are converted into a multiple integral over certain contours by means of appropriate 
functions. 
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A priori it is unclear how to choose these functions. They should be related to the 
functions appearing in the description of the thermodynamics of the spin chains. Yet, 
there are several mathematically rather different formulations of the thermodynamics 
using different types of auxiliary functions. In the study of the spin-^ chain [12, 13] 
only one of these formulations turned out to be compatible with the multiple integral 
representation. It is the formulation based on the quantum transfer matrix [31] and 
using only a finite number of auxiliary functions which satisfy a closed set of functional 
equations [25]. So far this is the least canonical formulation. No general scheme for it 
is known. Fortunately, the best understood case is just the case of the higher-spin XXX 
chains, which was worked out by one of the authors [30]. As we shall see below the 
auxiliary functions introduced in [30] are indeed most useful also in the framework of 
multiple integral representations. These functions can be efficiently calculated from 
a set of nonlinear coupled integral equations and allow for an accurate numerical 
description of the thermodynamics of the higher-spin chains [30]. Besides the auxiliary 
functions that satisfy nonlinear integral equations we shall introduce new functions, 
solving linear integral equations, which will finally allow us to rewrite the combinatorial 
sums representing the density matrix as a single multiple integral. 

We see this work as a feasibility study and therefore stick with the simplest higher- 
spin generalization of a finite-temperature multiple integral representation. Further 
generalizations to general higher spin, to the XXZ case or to include a disorder para- 
meter into the calculation are left for future studies. 

The paper is organized as follows. In section 2 we recall the construction of the 
Hamiltonian and the statistical operator by means of fusion of spin-^ transfer matrices. 
We also recall how to calculate the density matrix of a chain segment within the quantum 
transfer matrix approach. In section 3 we review the calculation of the thermodynamic 
quantities by means of nonlinear integral equations and present an alternative closed 
contour form of such equations. Section 4 contains our main result, which is a multiple 
integral formula for the inhomogeneous density matrix of a finite chain segment. In 
section 5 we present a factorized form of our formulae for the one -point functions. 
Finally, the zero temperature limit is sketched in section 6. The technical details of the 
derivation of the nonlinear integral equations and of the multiple integral formula have 
been separated from the main text and are summarized in three appendices. 

2 Hamiltonian and density matrix 
2.1 Hamiltonian 

The Hamiltonian of the integrable isotropic spin-1 chain on a lattice of 2L sites is 

H = J . t (St^-(St^) 2 ). (1) 

* n=-L+l 

Here implicit summation over a = x,y,z is understood, and periodic boundary condi- 
tions, 5" L = 5", are employed for the explicit sum over n. The 5™ act locally as standard 
spin-1 operators, and antiferromagnetic exchange, / > 0, is assumed throughout the 
paper. 
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The Hamiltonian (1) was first obtained in a more general anisotropic form in [34]. 
Shortly later it was constructed by means of the fusion procedure [26,27]. The ground 
state and the elementary excitations were studied in [33], and an algebraic Bethe ansatz 
and the thermodynamics within the TBA approach were obtained in [1]. 



2.2 Integrable structure 



The model can be constructed by means of the fusion procedure [26], starting from the 
fundamental spin-^ /^-matrix 



(\ 



\ 



b{X) c{X) 
c(X) b(k) 



\ 



V 



b{X) 



X + 2i' 



c(X) 



2i 
l + 2i 



(2) 



which we think of as an element of End(C 2 <8> C 2 ). It satisfies the Yang-Baxter equation 

R 1 V ] (X-H)R^ ] (X)R^ ] Oi) =R^\h)R^ ] (X)R^ ] (X-») . (3) 

As usual the R^ in this equation act on the j'th and kth factor of the triple tensor 

product C 2 C 2 ® C 2 as /ft 1 ' 1 ! and on the remaining factor trivially. PJ 1 ' 1 ! is normalized 
in such a way that 

R^(0)=PW, (4) 

where P W is the transposition of the two factors in C 2 ® C 2 . We say that PJ 1 ' 1 ' is regular. 
At the same time R^'^ = pW/?! 1 ' 1 ! satisfies the unitarity condition 



R^(X)R^ l] (-X)=I 4 , 



(5) 



with /„ denoting the n x n unit matrix. 

A further property of R^- l \ which is at the heart of the fusion procedure, is its 
degeneracy at two special points, 



lim = P ± 

X-^±2i 2b{k) 
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• (6) 



The P ± are the orthogonal projectors onto the singlet and triplet subspaces V^ s \ C 
C 2 <8> C 2 with standard bases 

b( ' 5) ={ti(U)^(?)-(?)^U))}' 

« (r) ={ (i) ® (i) . ^ ((i) ® (?) + (?) ® (5)) , (?) ® (?) } - a) 

Due to (3) and (6) we have the important relation 

P 23 R[ n ] &) r[ U ] + 2i ) P 2 + 3 = , (8) 
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meaning that R [ ^ l] (X)R [ ^ l] (X + 2i) leaves C 2 (g> VW invariant. 
Let us introduce U : C 2 <g> C 2 ->• C and 5 : C 2 (g> C 2 ->• C 3 , 



^=(0,^,-^,0), 5: 



/l 



V 



J_ _j_ 

y/2 V2 



(9) 



V 



which map the singlet and triplet subspaces of the tensor product of two spin- ^ repre- 
sentations onto C or C 3 , respectively. These matrices satisfy 



ss'=h, s's = p + , uu' = i, u'u = p- 

where the superscript t indicates the transposition of matrices. 
Using S we can define the fused /^-matrices 

(X) = S 23 R[ l 3 1] (X)r[ 1 2 1] (X + 2i) S 23 , 



rW(X) = 5 12 /?S 1 3 ' 11 (X-2i)4 1 3 11 (X)S[ 2 , 



[Ml, 



(10) 

(11a) 
(lib) 



R l 2 > 2 \l) = 5 12 534/?Si 1 k^-2i)/?S 1 3 1 k^)4yk^)4 1 3 11 (^ + 2i)5 f 345 f 12 (11c) 

acting on C 2 ® C 3 , C 3 ®C 2 , or C 3 ® C , respectively. Combining the Yang-Baxter 
equation (3) and equations (8), (10) it is easy to see that 

/?^' 2i2] - / /)/?s 2 3 si ' 2s3] (X)4 2 3 ^' 2s3] = 4 2 3 i2 ' 2s3] (.u)/?;^ 1 ' 2531 (x)/?^ 1 ' 2 ^ 1 (x- ju) , d2) 

where Sj = \, \ for j = 1,2,3. 

In particular, PJ 2 - 2 ! is a solution of the Yang-Baxter equation. With pPl denoting 
the transposition on C 3 (8) C 3 and pP> 2 l = P PlpJ 2 ' 2 ! it has the further properties 



R^(0)=pW, 
RW(X)RW(-X)=I 9 , 



(13a) 
(13b) 



i.e., PJ 2 ' 2 ' is regular and unitary. It follows with (13a) that P^ 2 ! generates the Hamilto- 
nian (1), 

i[ 2 , 2 ] 



H = iJ £ h n - hn , h n - hn = d x K-l n m x =o- 

n=-L+\ 



(14) 



2.3 Density matrix 

In [13] we have set up a formalism which enables us to calculate thermal correlation 
functions in integrable models with R- matrices fulfilling (13a). It is based on the so- 
called quantum transfer matrix [3 1 ] and its associated monodromy matrix which are 
directly related to the statistical operator. 

The Hamiltonian (1) preserves the total spin 

L 

S a = £ Sf. (15) 

j=-L+\ 
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Thus, the magnetization in z-direction is a thermodynamic quantity, and the statistical 
operator 

H-2hS z 

p L (T,h) = e-— (16) 

describes the spin chain (1) in thermal equilibrium at temperature T and magnetic 
field h. 

The statistical operator does not exist in the thermodynamic limit. Quantities that 
are better defined for the infinite chain are the free energy per lattice site and the density 
matrix of a finite chain segment. The free energy per lattice site is 

f(T,h) = -T hm lntr - L+1 "" XpL(r ' /i) . (17) 
l^oo 2L 

It determines the thermodynamics of the model [30] which will be briefly reviewed in 
section 3. The density matrix of a finite chain segment [l,m] is defined as 

Dm m] (T,h) = hm "' '"' 1 '^'''^ . (18) 
l ' 1 tr_ L+ i,..., L p L (r,/i) 

With D[! m ](T,h) we can calculate the expectation value of any local operator that 
acts trivially outside the finite segment [l,m]. In particular, D^ m ](T,h) allows us to 
calculate the static correlation functions inside [l,m]. 

For any integrable model, whose /^-matrix does not only satisfy the Yang-Baxter 
equation, but also the regularity and unitarity conditions (13), we can approximate the 
statistical operator Pl(T,K) of the 2L-site Hamiltonian using the monodromy matrix 
of an appropriately defined vertex model with 2L vertical lines (— L + 1 , . . . , L) and N 
alternating horizontal lines (1, . . . ,N with N even). This fact was exploited many times 
in the calculation of the bulk thermodynamic properties of integrable quantum chains, 
in particular, in case of the higher- spin integrable Heisenberg chains [30]. In [13] it was 
noticed that the same formalism is also useful for the calculation of thermal correlation 
functions. Following the general prescription in [13] we define 



T m {X) = e^/^gl (X - P/A0*J^(-P/tf " A.) • • • 

. . . R l2 ? ] (X - P/A04 2 ' 21 H (-P/W - X) , (19) 

where t\ indicates transposition with respect to the first space in a tensor product. This 
monodromy matrix is constructed in such a way that (see [13]) 

L 1 - 

*1...n{t%(0)...t[ 2] (0)}= [l- A _^(prA„-i )l ,-2^) + 0(^)J ' (20) 

Hence, setting P = U/T and 

p N ,L(T,h) = tr- L ^{T% (0) . . . T[ 2] (0)} , (21) 
we conclude, using (1), (14) and (20), that 

\imp N . L (T,h)=p L (T,h). (22) 
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We shall call this limit the Trotter limit. 
The transfer matrix 

t^{X)=tVjTf\X) (23) 

is commonly called the quantum transfer matrix. We shall recall below how it can be 
diagonalized by means of the algebraic Bethe ansatz [30]. Quite generally it has the 
remarkable property that the eigenvalue A PI (0) of largest modulus of f PI (0) (we call it 
the dominant eigenvalue) is real and non-degenerate and is separated by the rest of the 
spectrum by a gap [31, 32]. It can further be shown that 

f(T,h) = -T lim lim ^-L + ^xP N , L (T,h) = _ T ^ ^ ^ 

Thus, the dominant eigenvalue alone determines the bulk thermodynamic properties of 
the spin chain. 

Owing to the fact that R^ 2 ' 2 ^ satisfies the Yang-Baxter equation the transfer matrices 
?I 2 1(X) form a commutative family, 

[t®(X),t®( M )]=0. (25) 

It follows that the eigenvectors of t^(X) do not depend on X. Let |*Po) denote an 
eigenvector belonging to the dominant eigenvalue A' 2 ! (0). We shall call it the dominant 
eigenvector. It is unique up to normalization and is an eigenvector of t^(k) with 
eigenvalue A® (k) = k [2] (A.) I^Fo) / <^o |^o> • In [13] it was pointed out that such an 
eigenvector determines all static correlation functions at temperature T and magnetic 
field h. In particular, it determines the density matrix (18) of any finite segment [l,m], 

For technical reasons it is better to consider a slightly more general expression 
than the one under the limit, by allowing for mutually distinct spectral parameters 
j = l,...,m, instead of zero. Setting ^ = (^i,...,^ m ) we define 

D[ 2] rn = ^oirP^O^-^r^Ul^o) m , 

U ^ " <Vo>A[ 2 ^i)...A [2] (U ' 
the inhomogeneous density matrix at finite Trotter number. Then 

D [Lm] (T,h)=hm lim D^). (28) 

The expression (27) is our starting point for the derivation of the multiple integral 
representation in appendix C. 



2.4 Bethe Ansatz solution 



For the calculation of the free energy (24) and the inhomogeneous density matrix (27) 
we need to know in first place the dominant eigenvector |*Po) and the corresponding 
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transfer matrix eigenvalue A^(X). They can be obtained by means of the standard 
algebraic Bethe ansatz for the spin-^ generalized model (see e.g. chapter 12.1.6 of [11]), 
since, by the general reasoning of the fusion procedure [27], the quantum transfer 
matrix ?P1 (X) can be expressed in terms of a transfer matrix with spin-^ auxiliary space 
and its associated quantum determinant. 

For the temperature case at hand we define the staggered monodromy matrix with 
spin- j auxiliary space [30] by 

rj l] (x + i) = ^' t r [ i§ (i - m)R [ N-t(-m - a.) . . . 

. . . R l ^ ] (X - P/AOM 2 . 1 * (-W - *>) ■ (29) 

Then, interpreting this monodromy matrix as a 2 x 2 matrix in the auxiliary space a, we 
define 

t^(X)=trT^(X), det q T^(X) = U(T^(X-i)^T^(X + i))U t . (30) 

It follows from (11) that 

T^(X)=S(T^(X-i)^T^(X + i))S t . (31) 

Taking the trace and using (30) we conclude that 

t [2] (X) =?W(X-i> [11 (X + i)-det q rW(X), (32) 

Hence, since det q T^(X) commutes with T^(X) [27], every eigenstate of t^(X) is an 
eigenstate of t^(X) as well. 

The algebraic Bethe ansatz is based on the Yang-Baxter algebra relations 



which follow from (12) and (29). Representing T^(X) by the 2 x 2 matrix 



T^(X) 

and defining the pseudo vacuum 

10} = 



f A(X) B(X) 

MX) Dm 



(33) 



(34) 



L\l 



) <g> ( o 



we deduce from (29) that 

Cmj0>=0, A(X)\0)=a(X)\0), D(X)\0) = d(X)\0), 



(35) 



(36) 



where the pseudo vacuum eigenvalues a(X) and d(X) are explicit complex valued 
functions. Using the notation 



<\>±(X) = {X±iu) N ' 2 , u = - 



NT 



(37) 
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which proved to be useful in [30], we can express them as 

eV^-fl + i) e-*/ r 4> + (A,-i) ns . 

a{X)= ^_(X-3i) ' W= ^ + (^ + 3i) ■ (38) 

Given the Yang-Baxter algebra (33) and the pseudo vacuum eigenvalues (38) the 
eigenvectors and eigenvalues of can be obtained from general considerations 

(see e.g. chapter 12.1.6 of [11]). The dominant eigenstate \*¥o) of t^(X), in particular, 
can be represented as 

\Wo)=B(Xi)...B(X N )\0), (39) 

where the set of so-called Bethe roots {Xj}^ is a specific solution of the Bethe ansatz 
equations 

For the given set of Bethe roots {A,y}yL 1 we define the Q-function 

q (X) = f\(X-Xj). (41) 

7=1 

Then the eigenvalue of fW (X) corresponding to \Wo) is 

As for the eigenvalue of fPl (X) we conclude with (32) and equation (48) below that 

API (X) = AW (X - i)AW (X + i)-a(X + i)d{X - i) (43) 

This eigenvalue and the Bethe ansatz equations (40) are the main input for the calculation 
of the thermodynamics of the spin-1 chain. In order to perform the Trotter limit the 
eigenvalue must be represented by means of auxiliary functions satisfying a finite set 
of nonlinear integral equations. This was achieved in [30]. To the extend we need the 
results also for the calculation of the density matrix, they are reviewed in the following 
section. 

2.5 Simplified form of fused monodromy matrix 

Slight simplifications are possible for the form (31) of the fused monodromy matrix 
rPl and for the form (30) of the quantum determinant of jW. We include them here for 
later convenience. Setting r± = rW(^±i), 

I ' I (44) 




and using the Yang-Baxter equation and (10), we conclude that 

r [2] (^) = S(T + ®T~)S' =S(T~®T + )S'. (45) 
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Similarly, it follows that 



S(T + ® T~)U' = U(T~ <g> T + )S' = 
with the help of which we can represent e.g. as 



(46) 



(A-A+ V2A+B- B+B 

V2C-A+ C B++D A+ V2D B 
C C + V2C+D- D D+ 



) 



(47) 



anddetqrl 1 !^) as 



det q T^^)=D-A+ -B-C+ . 



(48) 



3 Thermodynamics 

In this section we consider the evaluation of the free energy per lattice site f(T,h) by 
means of nonlinear integral equations (NLIE). This gives us the opportunity to introduce 
certain auxiliary functions and integration contours that are also relevant for the multiple 
integral representation of the density matrix elements in the next section. Our starting 
point is the expression (24) for f(T,h) in terms of the dominant eigenvalue A PI (k) of 
the quantum transfer matrix together with the Bethe ansatz solution (40)-(43). In [30] 
the problem was solved within the more general context of the fusion hierarchy, and 
NLIE for the integrable isotropic spin chains of arbitrary spin were obtained. We believe 
that those NLIE are optimal in several respects for the calculation of the free energy. 
They are integral equations of convolution type formulated for a minimal number of 
functions on straight lines, and, for this reason, can be accurately solved numerically. 
Moreover, the low temperature asymptotics of the free energy can be extracted from 
these equations [30] . 

For the calculation of the free energy for spin 1 we will be dealing with three 
coupled NLIE for three functions b, b and y. We show the equations below in (58) 
and present an alternative derivation in appendix B. For finite Trotter number N the 
functions b, b and y can be expressed in terms of the Q-f unctions (41) and the functions 
(j)± introduced in (37) (see appendix A). This defines them as meromorphic functions 
in the entire complex plane, but is inappropriate for performing the Trotter limit. In the 
NLIE, on the other hand, the Trotter number appears only in the driving term and the 
Trotter limit is easily obtained. For a discussion of some of the subtleties related to the 
Trotter limit and the definition of useful auxiliary functions see [15]. 

If one is only interested in the free energy, it is sufficient to know the functions b, b 
and y close to the real axis (see (58), (62) below). For the calculation of more general 
physical quantities, however, as, for instance, the density matrix elements we are going 
to consider in the next section, we need to know b, b also close to straight lines parallel 
to the real axis, passing through ±2i. This is the reason why we reconsider and slightly 
extend the approach of [30]. 

The necessity of considering auxiliary functions in an extended strip around the 
real axis originates from the particular distribution of the Bethe roots that parameterize 
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^2j-l x 



-2i 



Im 

2i 



Figure 1: Schematic distribution of the upper and lower Bethe roots 'kij-x and %2j, 
respectively, in the strips § ± . 



the dominant state. Define the strips 

§ ± = {XGC|0<±ImX<2}. (49) 

Then the Bethe roots of the dominant state come in N/2 pairs (so-called two-strings) 
with one root in § + and the other one in S~. For large Trotter number they accumulate 
in the vicinity of ±i. We shall call the Bethe roots in § + the upper Bethe roots and the 
Bethe roots in S~ the lower Bethe roots. By convention the upper Bethe roots will be 
denoted A,2 7 - 1 and the lower Bethe roots X^j, where j = 1, . . . ,N/2 (see figure 1). 

Typical physical quantities at finite temperature can be written as sums over the 
Bethe roots of the dominant state. Such sums can be converted into contour integrals 
by means of appropriate auxiliary functions having their zeros at the Bethe roots. As 
compared to the spin- \ case the choice of the contours and auxiliary functions is more 
delicate for spin 1 . In particular, it seems that the auxiliary functions and integration 
contours have to be chosen separately in S + and S~. We shall consider the auxiliary 
functions 

w-ipib)' f(X) =f(^ (50) 

(see appendix A for the definitions of b and b in terms of ^-functions). As usually we 
also introduce the corresponding 'capital functions' 

$(X) = l + f(X), SM = l+fW. (51) 

They are meromorphic for finite Trotter number, and # has in § + exactly N /2 zeros 
located at the upper Bethe roots and only a single Af/2-fold pole at i — iu. Similarly, J 
has in S~ exactly N/2 zeros located at the lower Bethe roots and only a single A r /2-fold 
pole at — i + iu. 
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Figure 2: The contours C 1 * 1 in the strips encircle the upper and lower Bethe roots 
respectively and close at infinity. 



Using this information and the definitions of some additional useful auxiliary 
functions in terms of ^-functions (see appendix A) we obtain the following NLIE, 

% (A.) n ' 93 (X.) 

+ [ ^ K (X- M )\n!S(M)+ I ^K(X- M )\aS(M), A G 6+ , (52a) 
Je+ 27ti Je 2tci 

In H nuJW 

- f_ + ^K(X-v)\nm- j & J^.K(X-v)lnm, * € C" . (52b) 
Here we have introduced the kernel 



and the driving term 



m = 2h j)+ (X - 3iH_ (X - i)(|)_ (X + i)<j) + (X + 3i) 
Aj r + <|)_(X-3i)(|) + (X-i)<|) + (X + i)<|)_(X + 3i) 

^ y + + (54) 

Note that the only explicit Af-dependence is in the driving term df. And since this 
driving term has a simple Trotter limit, we conclude that the functions f and f have a 
Trotter limit as well. 
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Figure 3: For the regularization in the multiple integral representation the dashed lines 
show the relative positions of the contours 6*, S . 



The precise definition of the integration contours is slightly subtle. We illustrate it 
in figure 2. C + is a simple closed contour inside S + that encircles the upper Bethe roots. 
We may realize it as a large rectangle with upper edge slightly below 2i and lower edge 
slightly above the real axis. Similarly C~ must enclose the lower Bethe roots inside S~ 
and may also be taken as a large rectangle, now with lower edge slightly above — 2i and 
with upper edge slightly below the real axis. The bar in C means that the contours do 
not encircle the singularities originating from the kernel K(k). This prescription may 
be seen as an 'ie-regularization' of the kernel after the contour integral is decomposed 
into an integral over straight lines. Such type of regularization is needed because the 
kernel has poles at /u = A ± 2i which must not lie on the contours. Having in mind the 
multiple integral representation in the next section we prefer to realize it in the way 
sketched in figure 3, where S + — 2i inside C = C~ inside C + — 2i, and 'inside' means 
'infinitesimally narrower'. 

At first sight, (52a) and (52b) do not seem to be enough to fix the unknown functions, 
as the number of equations is smaller than that of the functions. In order to understand 
that they actually fix the functions f and f, let us simulate one step in the iterative 
scheme. Assume that an approximate estimation of f,f is already known. Then an, an 
are determined from f , f by 

f(X+2i) ~~ f M . , „_ _ s f(X-2i) ~~ f ^ . . „+ 

an (A,) = - — =— , for A G C , a n (A) = - — -frr- , for A e C . (55) 



Note that 53 (A) /#( A) and 03 (A) /#( A) are equal to 1 + a n (A) and 1 + a n (A) , respec- 
tively. They are thus determined by given f , f. Substituting them into the rhs of (52a) and 
(52b) (and an, an into the lhs), we obtain the next-step approximation to f,f. Therefore 
equations (52) consistently fix f and f. The other functions are then determined from 
them. 

Suppose that we have evaluated the auxiliary functions through (52). Then, for 
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\lmX\ < 1, the largest eigenvalue A^(X) is obtained as 

lnA^(X)=lnA^ ] (X)+ f ^K(X-n~ 3i)ln<8(/i) + f ^(X- ^ + i)ln^) , 
Je- 2tci Je+ 27ii 

lnA' 2 i(X) = 2/i + ln^^- 2i ^^ + 2i ; ^ 2h + 4J ^. (56) 

oVy r ^_(X-2i)^+(X+2i) r rl 2 +4 v ' 

The NLIE (52) are actually only one of many possible choices. We choose this one 
as we think that it has an advantage compared to others in the following sense. Although 
the equations themselves are literally correct, the integrations over contours suffer from 
poor numerical accuracy, especially in the low temperature regime. Therefore it is 
better to rewrite them in the form obtained in [30], where the integrations are defined 
on the straight lines. We will show in appendix B that (52) can be transformed into (58) 
below with the help of additional algebraic relations among the auxiliary functions. In 
the same appendix B we also provide subsidiary equations that determine the functions 
f, f on straight lines close to the real axis, which amounts to knowing b and b on straight 
lines close to ±2i (see (50)). 

Unlike in (52) we need to deal with b(X), b(X) and y(X), if we choose straight lines 
as integration contours. For convenience we introduce the shifted functions 

b E (X) = b(X-ie), b E (X) = b(X + ie), (57) 

and similar capital functions. Then the desired NLIE read 

logy(X)\ ( \ / logy (A,) 

logb E (X) = A b (X) \ +%*{ log«B e (X) I , (58) 



logb £ (X)J \A F (A,)/ \log $ B e (X) / 
where (%*g)i denotes the matrix convolution £y f™ m d]uXij(X — n)gj(jj), and 

A b (X) = -^+d(u,X-ie), A B {X) = j+d(u,X + iz), (59a) 
,/ an N f°° „ in sinb.uk J 71 

rf („ A ) = -y_d^— _> - fJ ^j- r <5%) 

The integration constants (±h/T) are fixed by comparing the asymptotic values of both 
sides of (58) for \X\ — > °°. The kernel matrix is given by 



%(X) = 



( X(X+ie) X(X-ie) \ 

X(X-ie) 5(X) -J(X + 2i(l-s)) 

v X(X + ie) -J(X-2i(l-s)) ?(X) j 



(60) 



where 

dk e -W- ikX 



1 f 
% ^ ~ 4cosh;cV2 ' ~ J -co 271 2 cosh £ 



(61) 
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The free energy then follows from (24) noticing that the dominant eigenvalue can be 
represented by integration over straight lines as 

InA^mA^)-^ 

/oo /»co 
d i u3<:(X- i u + ie)ln5S e ( i u)+ / dpX(X- p-ie)^^) . (62) 
-oo J — oo 

As the actual transformation from (52) to (58) is involved, we defer the details to 
appendix B. 

4 The multiple integral representation 

In this section we present the main result of this work, which is a multiple integral 
formula for the matrix elements d 2 ^ 1 ' 'r™(^), «/, P& = — ,0, +, of the inhomogeneous 
density matrix (27). Our formula generalizes the result of [21] to finite temperature 
and magnetic field and the result of [14] to spin 1. The details of the derivation can be 
found in appendix C. 

For any two sequences (a) = (oc„)™ =1 and (|3) = (Pn)™=i of upper and lower matrix 
indices we shall obtain a different multiple integral. Let us introduce the notation n a (x), 
o = — ,0, +, (x) = (a), (|3), for the number of as in the sequence (x), e.g. «o(P) is the 
number of zeros in (P). Then 

n + (a) + n (a) +n_(a) = m, (63a) 
n+(p)+n (P) + n-(P)=m, (63b) 
n+(p)-n_(p)-n+(a)+n_(a) =0. (63c) 

Here the last equation is equivalent to 2n + (a) + no(a) = 2n+(P) + «o(P)^- 

The dependence of the multiple integral on the indices ay, P^ enters through a 
sequence (z) = (z n )lZi encoding the positions of — ,0,+ in (a) and (P). For the 
construction of (z) we order the density matrix indices as CC m , . . . ,CCi, Pi , . . . , P m and 
inspect them starting from the left. If a m = — we do nothing, if a m = we define 
Z\ = m, and if a m = + we define z\ = Zi = m. We continue this procedure with a m _i 
and so on. When we have reached (Xi we have defined 

p = 2n+(a)+n (a) (64) 

elements of the sequence (z) in this way. If Pi = — we define z p +i = z p +i = 1, if 
Pi = we define z p +\ = 1, and if Pi = + we do nothing. We continue the same 
way with P2, P3 etc. until we end at P, M . The sequence (z) thus constructed has 
2n+(a) + no(a) +«o(P) +2n-(P) = 2m elements, and the pair (z), p is in one-to- 
one correspondence with the sequences (a) and (P). As an example let us consider 
(a) = (+, - 0), p = (0,0, 0). Then zi = 3, zi = z 3 = U = 1, z 5 = 2, z 6 = 3, p = 3. 



1 Using (47) this translates into the fact that number of plus signs in the sequences of upper and lower 
indices of the matrices the density matrix element (C.25) is composed of, must be the same. 
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Two types of functions occur under the multiple integral. One type is explicit and 
has its origin in the Yang-Baxter algebra. The functions 

m i—\ m 

F e (X) = Yl(X-$ k -i)H(X-$ k -3i) n (*-k + i), (65a) 

fc=l k=\ k=e+i 

m i—\ m 

F ( iX) = Yl(X-^ + i)H(X-^ + 3i) f[ (^-^-i) (65b) 

k=i k=\ k=e+i 

belong to this type. We think of them as 'fused wave functions'. 

The other type is related to the task of rewriting sums over Bethe roots as integrals 
over closed contours (see appendix C). These functions may be defined as solutions 
of linear integral equations over closed contours. We have two pairs of such functions. 
The first one is defined by 

G+(X,(,)=K(X-%-3i)-K(X-i-i) 

G (X, ^)=K(X-^ + i)-K(X-^ + 3i) 

where X G S + for G + and X G C~ for G . The second pair of auxiliary functions needed 
in the definition of the multiple integral is 

S+(X,^ = - e (X-^-5i)- e (X-^-i)-^{K(X-^-3i)+K(X-^-i)) 
S-(X^) = - e (X-^-i)- e (X-^ + 3i)-^(K(X-^ + i)+K(X-^ + 3i)) 

-^wf^-^i^m^^ (67b) 

where, similar to the above case, X G S + for S + and X G 6~ for S~ and where we have 
introduced the 'bare energy function' 

The functions G ± and 5 ,± enter the multiple integral through the determinant of a 
matrix with elements &[ p ) defined by 

W jG+{^k) j=h...,P 
J > 2k ~ l G-fco^y j = p+l,...,2m, 
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(p) 
j,2k 



iS + ((Oj,^ k ) j=l,...,p 
iS~(®j,Z>k) j = p+l,...,2m. 



(69b) 



Using all of the above defined notation we can write the non-vanishing matrix 
elements of the inhomogeneous spin- 1 density matrix as 



2— m— n + (a)— n_ (p) 



ni<;<,< m fe-^) 2 [fe-y 2 +4] 



'P f dCO, , J [ r d(0, , , 



det 2m 0g 



- n.l<j<k<2m(®j-®k-2i) ' 



(70) 



This formula is the main result of our work. It represents the inhomogeneous density 
matrix of the integrable spin-1 chain as a single multiple integral. All dependence on 
the Trotter number has been absorbed into the auxiliary functions G ± and S ± . Therefore 
the Trotter limit is trivial in this formulation. 

Note that it is also easy to perform the homogenous limit. In complete analogy with 
the spin- \ case [14,24] we obtain 



Di 



«l,...,0C m ^j. ^ _ 2~ m n +( a )-' l -(P) 

2m ' d(0; 



" p rfto 



n 



e 27ti 



F Zj ((Oj) 



n 



-Fzj(®j) 



det 2m sg 



n 



l<j<k<2m 



for the physical density matrix. Here we introduced the notation 



(co ; --cOi-2i) 



(71) 



s (p) 



"7,2* 



G + (<0/^)k=o j = l,...,p 



(*-l)! lG-(co 7 ^)|^o y = p+l,...,2m, 



la*- 1 



I 5+(COy,^)|^o 7 = l,...,/> 

(*-l)! \5-(co y ^)|^ = o j = p+l,...,2m. 



(72a) 



(72b) 



5 One-point functions in factorized form 

In this section we have a closer look at the one -point functions which are the most 
elementary correlation functions. Using the general multiple integral formula (70) we 
can write the non-zero one-point functions as 



4 Je 2ni Je 2ni 



d(0 2 (COi 


-\- 


-i)(co 2 -£-i) 


G+(C0!^) 




2ni 


(Oi 


- C0 2 - 2i 




5 + (co 2 ,^) 


du) 2 (coi 




-i)(co 2 -£ + i) 






2%i 


(Oi 


- oo 2 - 2i 


G-((02,y 




doo 2 (coi 


-s- 


H)(ro 2 -^ + i) 


G-(C0!^) 


5-(co!,^ 


2%i 


COi 


- C0 2 - 2i 


G-(C02,^) 


5-(co 2 ,^) 



(73) 
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This is the double integral form of the one-point functions. Judging from our experience 
with the spin-j case [3,7] and with the spin-1 ground state correlation functions [21] 
we expect these integrals to factorize into sums over products of single integrals. 

This is indeed the case. For n = 0, 1 we introduce the following functions repre- 
sented by single integrals, 



" vs; 7e+27tiQS(X) 7e-27iiQ5(X) v ' s; ' 



(74a) 
(74b) 



/e+27ii5S(X) ' h- 2ni<S(k) 

Then, using tricks similar to those employed in [3], we obtain the 'magnetization' 

= a © (75) 
and the 'probability for measuring zero for the z-component of the spin' , 

Go© <Tl©-2i 



(76) 



in factorized form. They determine all one-point functions because of the relation 

DX®+ti a ®+DZ($) = l. (77) 

Note that Go and also the whole determinant in (76) must vanish for symmetry reasons 
if the magnetic field is switched off. 

6 The zero temperature limit at vanishing magnetic field 

All dependence on temperature of the multiple integral formula (70) is hidden in the 
functions G 1 * 1 and S ± . We obtain the ground state result for vanishing magnetic field 
by replacing these functions by their corresponding limits which have to be calculated 
from (66), (67). 

The temperature enters these equations through the functions b, b, f, f and y. How 
do they behave in the limit? We first look at the nonlinear integral equations (58). As 
T — > for h = 0, the driving terms in the equations for b and b both go to minus infinity 
pointwise. It follows that b, b — > on lines slightly below or slightly above the real 
axis. From the equation for y we conclude that y — > 1 close to the real axis. Then by 
equation (B.16) also y(X±i) — ^ 1 for X close to the real axis, and, using (B.15), we find 
that f,f ->-0. Thus, 



(78) 
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for X slightly below or above the real axis. The behaviour of these functions close to 
the lower edge of C~ and close to the upper edge of C + then follows from (A.9b), 

UmUm ^-*> 1. limlim ffl+ 2 » = l (79) 
r^o+/w()<B(X-2i) 2 r^o+^o«B(X + 2i) 2 

for X slightly above or below the real axis. 

Inserting (78) and (79) into (66) and (67) and using that F(^) — > 2 we obtain a set 
of linear integral equations of convolution type that can be solved by means of Fourier 
transformation. Some care is required with the relative location of the contours, though. 
Referring to the notation 

Z ++ (X&) = lim lim limZ+a + 2i-i£,a, (80a) 

e-s-O+r-s-O+fc-i-O 

Z + - (!,£,)= lim lim limZ+a + is,^) , (80b) 

Z- + (X&)= lim lim limZ-a-i£,B, (80c) 

e-s-O+r-s-O+A-s-0 

Z~ (A,,E)= lim lim limZ-(A,-2i + ie,E), (80d) 

e-s-O+r-s-O+A-s-0 

where Z = G or Z = S, we obtain the following results 

G ++ (U) = G-(U)=0, (81a) 

^» = -^« = 2d^5))' (81b) 

^^ =s ~^=^m=wy (81c) 

J- + (U) = 7t( ^" 2 & \\ , ^(U)= ^~* + ff - (81d) 
v 4ch(f (*,-£)) ' v ' v 4ch(f (*,-£)) 

As a first consistency test we may insert these results into our formulae (74) for the 
one-point functions. We obtain (Jo = and 8o = —2. Then (75), (76) and (77) imply 
that £>+(0) = Dq(0) = Dl (0) = 1/3 as it must be from symmetry considerations. This 
is, of course, in agreement with [21]. 

Still, it is not obvious how to relate, in general, the limit of our multiple integral to 
the multiple integral derived there directly for the ground state at vanishing field. Here 
we consider only the case of the one-point functions and defer any further discussion to 
future work. We have to calculate the limits of G + and S + in the lower strip S~ and 
the limits of G and S~ in the upper strip S + on lines close to the real axis and close to 
±2i. These lines must be chosen in such a way that all poles of the kernels in (66), (67) 
are located outside the integration contours. Keeping this in mind we define 

Z ++ a,n= lim lim UmZ + a-i8,^), (82a) 
e^0+r^0+/i->0 

Z,(X,l)= lim lim limZ + (A,-2i + ie,E) , (82b) 
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Z-+Ck,l)= lim lim limZ~(A, + 2i-ie,E), 

e-s-O+T-s-O+fc-K) 

Z(X,l)= lim lim limZ-(A, + ie,E), 

e-s-O+T-s-O+fc-s-O 

for Z = G and Z = S. Inserting (81) into (66), (67) we obtain 

G_ + (U) = G + _(U) = 0, 
S_ + (U)=S+-(U) = 0, 
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2ch(f(A,-y)' 



7t(X-^-2i) -7C(X-^ + 2i) 

5 - (A " §)_ 4ch(l(A,-5))' ++( ' g "4ch(f(^))- 
Inserting (81) and (83) into (73), in turn, we arrive at 

d+(o) =d:(o) = 



(82c) 
(82d) 

(83a) 
(83b) 

(83c) 
(83d) 



dxi 



4 j-oo ch(Ttxi) j-co ch(;tX2 



dx? 



(xi + j)(x 2 -h) (xi + j)(x 2 - 



"4 

71 

2i-, 



dxi 



ch(TtXi) 

x 2 + 



dx? 



dx 



, ch(7tX2) 

1 



x\ — X2 + i0 xi — X2 — iO 

(Xi + ^)(X2-|) {xi-^){x2+{] 



x\ -x 2 + i0 



x\ — X2 — 2i 



ch 2 (;tx) 3 ' 



(84) 



to be compared with (4.9) and (4.13) of Kitanine [21]. Similarly, (4.12) of [21] for 
D[J(0) is reproduced as well, 



D § ( o)4£ 

_7t r 



dxi 



dx? 



ch(Ttxi) j-oo ch(;tX2) 



(x! + ^)(x 2 + ^) (xi + i)(x 2 + i) 



x\ — X2 — iO 



xi -x 2 +i0 



dx 



-2x 2 
ch 2 (7tx) 



1 



(85) 



7 Conclusion 

We have managed to represent the inhomogeneous density matrix of the integrable 
isotropic spin-1 chain as a single multiple integral (70). Our formula admits of the 
Trotter limit, the homogeneous limit and the zero temperature and zero magnetic field 
limit, where it reproduces the known values of the one-point functions. The main 
difficulty in the derivation of (70) was not in the algebraic part, which can be treated in 
a similar way as in the ground state case, but in the analytic part. For finite temperature 
we can not work with root density functions. Instead, the integrals are obtained by 
replacing sums over Bethe roots by integrals over closed contours encircling the Bethe 
roots. In the spin-1 case the Bethe roots for the dominant state of the quantum transfer 
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matrix come in widely separated pairs, so-called two-strings. In the Trotter limit they 
cluster close to ±i. Therefore, in order to avoid unwanted extra-terms, we were forced 
to introduce closed contours consisting of two separated loops, which brought about 
a considerable amount of technical complexity into the derivation as compared to the 
spin- j case [12] (see appendix C). 

We believe that our result can be generalized to the critical anisotropic case, as it 
was done for the ground state at vanishing magnetic field in [10], and to arbitrary higher 
spins. Of particular interest for our own research will be the question if the correlation 
functions of the integrable higher spin chains factorize. We have obtained a first hint in 
this direction: we saw in section 5 that the integrals for the one-point functions factorize. 
This is still not what was called factorization of correlation functions in [4] and what was 
recently proved to hold for the spin-^ XXZ chain, namely, that all correlation functions 
(of a suitably regularized model) can be expressed in terms of a small number of special 
short-range correlations functions constituting the 'physical part' of the problem (for 
the physical part of the XXZ spin-^ correlation functions see [2]). Showing this for 
the higher- spin chains of fusion type as well will be a challenging project for future 
research. 
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Appendix A: Auxiliary functions for spin 1 

As long as the Trotter number is finite the transfer matrix eigenvalues AM (X) and 
A^ (X) as well as all the auxiliary functions used in this work can be expressed in 
terms of the ^-functions (41) and the functions ty± defined in (37). In this appendix we 
collect the corresponding formula and also some of the relations between the auxiliary 
functions. The presentation largely follows [30]. 

It is sometimes more convenient to deal with polynomials rather than with rational 
functions. For this reason a different normalization of the elementary /^-matrix was 
used in [30]. This leads to differently normalized transfer matrix eigenvalues. In order 
to simplify the comparison with [30] we define the functions 

Ai(X) = <b-(X-3i)§+{X + 3i)A [[] (X) , (A.la) 
A 2 (X) = <|>_ (X - 4i)<|> + {X + 2i)<|>_ {X - 2i)<|> + {X + 4i) A^ (X) . (A. lb) 

Then, following [30], we introduce 

h (X) = e" * <>_ - 4i)4>+ (X - 2i)<>_ (X - 2i)<|> + (X) , (A.2a) 

X 2 (X) = §-(k- 2i)<> + (X + 2i) g^j|M^±|j , (A.2b) 

X 3 (A) =e^_ (X>+ (X + 2i)4>_ (X + 2i)(|> + (X + 4i) q %P^ . (A.2c) 
It follows that 

A 2 (X) = X l (X)+X 2 (X)+X 3 (X). (A3) 
The basic auxiliary functions for spin 1 are 

X l (X + i)+X 2 (X + i) - x 2 (X-i)+X 3 (X-i) 
b(X) = X3IX+Y) ' b(X) = XAX^i) ' (A " 4) 

with corresponding capital functions 



05(A) = l + b(X), ( B(X) = l + b(X). (A.5) 

In [30] the nonlinear integral equations (58) were derived from a set of functional 
equations satisfied by the functions b, b, 23, 53 together with 

In appendix B we present an alternative derivation starting from the integral equations 
(52) and combining them with some of the algebraic relations exposed below. 
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In the derivation of the multiple integral representation for the density matrix 
elements we further encounter the functions 

n W = ^T> 2l(A) = l + a(A) (A.7b) 
a(A) 

familiar from the spin-j case. We find it also convenient to give a separate name to the 
functions with shifted arguments, 

a n (A) = a(A + 2i), a n (A) = a(A-2i). (A. 8) 

The following relations among the functions are needed at several instances in this 
work. They follow directly from the above definitions, 

S^ = 2l(A + 2i), f^=2l(A-2i), (A.9a) 
TO) V ^ TO) 

%(X + i)W(X-i) = l+y(k), (A.9b) 

<B(A-i) Agflj _ 1 

2l(A + i)2l(A-i) ~ AW(A + i)AW(A-i) ~ l+y-^X)' 

b(k-i)b(k + i) = \+y(X) = Y(k), (A.9d) 



(A.9c) 



TO »(A,) y(A,+i) 

TO = ^)F(^T) a(X) ' (A " 9e) 

an (A) \ 03(a) /y(A + i) 

/5(A) \ TO) 

^TO = (l(A) ( )= ^1^)- ( g) 



Appendix B: NLIE with straight contour integrations 

In this appendix we will show the steps that are necessary for transforming (52) into (58). 
We also present subsidiary equations which can be used for the numerical calculation 
of some of the auxiliary functions on lines away from the real axis. 

First note that numerical calculations with fixed Trotter number N suggest that 



|b(A)|, |f(A)|<l forImA = £, \b(k)\, |f(A)| < 1 forImA=-£, (B.l) 

\b(X)\, |f(A)| > 1 forImA = 2-£, \b(7i)\, |f(A)| » 1 forImA=-2 + £ 
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in the low temperature regime. Therefore we rewrite, for example, 



for X, located inside a narrow strip So including the real axis. To emphasize the relative 
location of X and (u, we write the last integral as 

-f ^K(X- M -2i)ln^4. (B.3) 

Jlmtelm/j 27C1 b{fi) 

We keep our assumption that X G So for a while. Thanks to (A.9f), (A.9g) and a 
similar transformation applied to the integrands, (52) is represented as 

lny(X - i) = -d f (X) + In® (A,) 

JImX>Im^27Il D(/lJ JlmX<lm/j 2711 b(/l) 

+ / ^T(A-/i+ie)lnSGu-ie)- / ^^(X-/i-ie)ln^(/i+ie), (B.4a) 



lny(A + i) = d/(A)+ln<8(A) 

-/ ^-,-2i)lnp+/ ^-, + 2i)ln5M 

-/ p T K(X-v + i£)ln$(n-i£)+ f ^^(X-/i-ie)ln^(/i+ie) . (B.4b) 
J-c»27Cl J-oo27Cl 

The integrands in the last two terms in (B.4a) and (B.4b) become proportional to the 
logarithm of 3(X)/$(X) in £ — > limit. Since such ratio does not appear in (58), we 
would like to replace it using (A.9e). For this purpose, we first note a contour integral 
representation for lna(A), 

lna(A) = -d f (X)-f e+ ^K(k-n)\n^)-J e _ ^(X-^ln^u) . (B.5) 

Again we rewrite this using integration on straight lines and substitute the result into 
(A.9e). It is then immediately clear that 



JImX>Im^27Cl b{fJ) JImX<Im^27Cl b(fJ) 

+ [ ^T(A-/i + ie)lnS(/i-ie)- [ ^U(A-/i-ie)lnSGu + ie) . (B.6) 
J-oo 2711 J-oo 27C1 
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To proceed further, it is convenient to consider equations in Fourier space. For a 
smooth function f(k) we define 

We also introduce shifted functions 

b e (A,) = b(A,-ie), b e (X) = b(X + is) , (B.8a) 

f e (X) = f(X + is), f e (A,) = f(A,-ie), (B.8b) 

and similarly for the capital functions. 

First we take the Fourier transformation of (A.9d) for X real. This leads to a direct 
relation between dlb e (&) and dlb e (&), 

e-^dlb e (it) = -e- {2 -^ k d\b e (k)+e- k d\Y(k)-iNe- k smhuk. (B.9) 

Similarly, take the Fourier transformation of (B.6) and delete dlb e (fc) by means of (B.9). 
Then 

e-*<&&(*) -e^dlUk) = y+ l Ki{k) {-dlAf(k) + e £k (K yh (k) - \)dm e (k) 
-e- £k (K yh (k) - l)m E (k)-^(K yb (k)+ e - 2k K y - b (k))d\b e (k) 
+ (t- k K yh (k) + t k - e"*)<fly(*)) , (B.10) 

dlA/(it) =We- 2 W-*sinh M A:, 

K 2 (k)=e-W, K yh {k)=z- m+k \ K yh (k) = e- 2 ^- k K 

Finally take the Fourier transformation of the logarithmic derivatives of both sides 
of (B.4a) and (B.4b). Note that d\$ £ (k) and d\$ e (k) only appear in the combination 
e- sk dl$ E (k) - e ek d\$ e (k). Therefore, by substituting (B.9) and (B.10), one obtains two 
equations containing dlb e , dly, dl5S e ,dl5S e and diy. They can be solved for dlb e and 
dly in terms of dl*B e , dl*B e and dly, yielding 

iMk) _ -n'-"*** 4- ^§ - «E§ e *-e,» + , (B. 1 la) 

v 7 2cosh£ l + e 2 W l+e 2 W 2cosh£ 

*W = ^^#t£^ffi. (B , lb) 

/v ; 2coshfc 

If dlb E (&) is eliminated from (B.lla) by means of (B.9), an equation for dlb E (fc) is 
obtained, 

3S. W = -iN e 4^ + - ^|e-^-' + fa t .. (B. 12) 

v 7 2cosh£ l+e 2 W l+e 2 W 2cosh£ 

Applying the inverse Fourier transformation and integrating once, we successfully 
recover the NLIE (58) with straight integration contours. 
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To evaluate physical quantities beyond A PI, we also need NLIE (denned with 
straight integration contours) for f E and f E . This can be understood as follows. The 
eigenvalues of physical quantities are parameterized by BAE roots. Thus, they can 
be naturally represented by loop integrals involving *8(X),*B(X),$(X) or $(X). We 
consider, for example, 

where P(X,n) is some function. This integral can be represented as 

1= -r d ,q^ + r d /^--*+ i£ > . (B.14) 

We therefore need to evaluate f e (X) and f e (A,), when we adopt straight lines near the 
real axis as integration contours. 

Indeed, it is not difficult to derive the following expressions for f e (X) and f e (X,), 

/oo 
^3C fe (A,-//)ln?8 e 0i) 
-oo 

/OO /''OO 
d/ift H (A,-/i)ln» e (/i)+/ d/i3C /3? (X-/i)lny_Cu), (B.15a) 
-oo J — oo 

/oo 
d/iDC^-^ln^Gu) 
-oo 

/OO /"'OO 
d f tX lA (X- f t)lnX E ( M )+ dnX fy {X-n)\nY + (n). (B.15b) 
-oo J — oo 



+ 



The integration kernel 3C a( >, is the corresponding component in (60), except for Xf y and 
X fy , defined explicitly by X fy (X) = -X (X - i ( 1 - e) ) and X fy {X) = - X (X + i ( 1 - e) ) . 

The functions Y±(X) denote shifted F-functions, Y±(X) = Y(X±i). Unfortunately, 
they can not be determined from (58), because of the singularity of the kernel function. 
We thus need subsidiary equations, 

lny + (X) = ln^ e (X + ie)+ [ d/i£j*(A,-/i + i)lni8 e (/i) 

J — oo 

/oo 
d M X yi (X-n + i)lnm e (n), (B.16a) 
-oo 

/oo 
d^^j*(A,-/i-i)ln?8 e (/i) 
-oo 

/oo 
d^^(A,-/i-i)ln® e 0i). (B.16b) 
-oo 

The functions b(X) and b(X) are analytic in a narrow strip including the real axis. For 
this reason we can use (58) to estimate the first terms in the rhs of (B.16). Thus, (B. 15) 
and (B.16), together with (58), fix f(X) and f(X) through integrals defined on straight 
contours. 
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Appendix C: Derivation of the multiple integral representation 

In this appendix we derive the multiple integral representation of section 4. Our strategy 
is to use as much as possible the results obtained in [12] for the spin-1/2 case. 

C.l Results for spin-1/2 auxiliary space 
C.l.l Spin projection conserving basis 

The monodromy matrix preserves the pseudo spin projection 

if = £(-i)^., [ri 1] (X),i^+^]=o. (c.i) 

It follows that 

T^)... T m;:(w (c.2) 

Since the dominant state |*Fo) = . ..B(Xn)\0) has pseudo spin projection zero, 
if |¥ ) =0, we conclude that the matrix elements (^ol^ 1 ^ 1 . . . rWjJ^)^) all 
vanish, unless £" =1 (aj — J3_/) = 0. 

This means that we must have the same number of plus signs in the sequences (a,-) 
and (Pfc) of upper and lower indices. Let us introduce a basis on the space of local 
operators which is adapted to this fact. It is convenient to label the states in this basis 
by the positions of the plus signs in (ay) and minus signs in (p^). For x = (xi , . . . ,x„) 
with Xj G Z n = {1, . . . ,n} and {xi , . . . ,x p }, {x p+ \ ,x n } two sets of mutually distinct 
numbers, let 

M x ) = °x n ■ ■■°x p+1 Gi ■ ■ *t<5- p ...a~. (C.3) 

Then 

+ if j G {xi,...,x p } 
— else 

b p (x)=e^...e n ^ with v ' (C.4) 



+ ifj<£{x p +i,...,x n } 
else. 



Clearly 

B n = {b p (x)\n > x\ > ■ ■ ■ > x p > 1 < x p+ i <■■■ <x n <n;p = 0,...,n} (C.5) 
is a basis of the r\ z = subspace of the space of local operators acting on (C 2 ) <8 "\ 

C.1.2 Combinatorial formula for density matrix at finite Trotter number 

Referring to the notation of the previous subsection we now fix an even n = 2m and a 
vector x that specifies a basis element in B2 m - We further define C, = . . . , C,2m) and 

n[% In _ (^o|tr{r[ 1 ](C 1 )0---0rW(UH(x)}|^o) 

° - (*o|^o)A[ 1 ](Ci)---AW(U) ' { } 
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Density matrix elements of this form were considered in [12], where a multiple integral 
representation for the spin- 1/2 XXZ chain at finite temperature was derived. Most of 
that calculation, up to the very last step, was purely algebraic and entirely based on the 
commutation relations between the elements of the monodromy matrix. This means 
it only depended on the structure of the /^-matrix and, hence, can be taken over to the 
present case. 

For this purpose let us first of all recall some of the notation of [12], but in a form 
already adapted to the rational limit. Let 



2m 



Fj l] (X)=ll(X-^-2i) n -CO, j=h-,P, 

k=xj+l 



(C.7a) 



2m 



Ff(X)=l\(X-^ + 2i) H (*,-£*), j = p+l,...,2m, (C.7b) 

k=xj+l 



k=\ 



and define a set of functions Wj(Q, j = 1, . . . ,N, as the solutions of the linear system 

a'(Xj)wj(Q = e(C - Xj)a(Q - - Q + f K(kj - h)MQ , (C.8) 

it=i 

where a is the auxiliary function defined in (A.7a). 

Then, from equation (63) of [12], we have the following combinatorial expression 



DW(x|Q= £ 

({ E +},{e-})G P2 (Z 2ra ) 



N 



2m 



+ .-4+ =1 (m,{S-})e P2 (Z 2m ) 

1 2m -" card{8-}=n 

sign(/?)det[-w £ + 



sign(Pe) 



Rg6 



IIl</<:*<2m(C* - 0)( ro 7 ~ ro * ~ 2i ) 



'*7 



2m— n 



n 



2m— n 



n -*S(kt 



7=1 



7=1 ' 



• (C.9) 



For the sums we have adopted the notation from [12]. 2,2m = {!,-• -,2m}, and pii^im) 
is the set of all partitions of Z2 m into ordered pairs of disjoint subsets. E.g. the first sum 
is over all pairs ({£+}, with U = Z 2m and {£+} n {er} = 0. Moreover, 
n = card{s~} by definition. We enumerate the elements in the sets {e^} and {8*} 
in such a way that zj < and 8^ < 8^ if j < k. Then for every ({£+}, and 
({8 + }, {8~}) the permutations P, Q G & 2m under the sum are fixed by 



Pj 



j= 



Qj 



5 7 
8 U 



j=\,...,n 
j = n + 1 , . . . , 2m . 



y e j-n j = n + l,...,2m 
Note that the Bethe equations a(kj) = — 1 were used in the derivation of (C.9) 



(C.IO) 



29 



C.2 Fusion for density matrix elements 
C.2.1 The narrow contour 

We shall employ equation (C.9) in the derivation of a multiple integral representation for 
the density matrix elements of the spin-1 chain. We begin by fixing real inhomogeneity 
parameters j = 1 , . . . , m, and £ > 0. For ;=l,...,mwe choose 8/ € { 1 + £, — 1 — e} 
arbitrarily and define 

^/ I C; • 18;. o; C/ 18/. (C.ll) 

Using the fusion formulae (31) and (45) we can express the inhomogeneous density 
matrix (27) of the spin-1 chain as 



Z)P](5) = lim 



s®»0Fo|r [1] (Ci) ® ■ ■ ■ ® rw (U)|Vo}(sT m 



(^o|^o)A[ 2 ](^i)---A[ 2 ](U 

: __ (Vol tr{rw (d) ® - ® rw (U)^(x)}| v ) 



= V S® m fc p (x)(S'f m lim 



= £ 5® m ^(x)(5 f )® m n 
xes 2 m L/=i 



(V |V )AP]£ 1 )...AP](J; m ) 
AW&-i)ANfo + i) 



AP](^.) 



limD [11 (x|Q. (C.12) 



Let us denote the coefficient under the sum 



n 

•7=1 



AW^-QAW^ + i) 
APl(^) 



lim£»W(x|Q. 

e->0+ 



(C.13) 



Inserting equations (C.9) and (A.9c) on the right hand side we obtain 



D[ 2 l(x|Q= lim £ £ 

fc '" ' )+ ({e+},{ E -})G P2 (Z 2m ) /<+,.••■<«+ =1 ({5+},{8-})ep 2 (Z 2m ) 

card{8 }=n 



III 

1 2m— n 

sign(^)det[-w/ + (Cs+)] 



sign(Pg) 



Jte6» IIl<j<*<2m(C* - - <B* - 2i) 



2m— n 



n 

7=1 7 7 



2m— n 



n -^(Ht 

7=1 7 7 



ii 

7 = 1 7 



n 

7=1 7 



(C.14) 
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Rc 




Figure 4: The contours P* 1 only contain Bethe roots {Xk} = {^-2/} U {X,2 7 _i} and 
the corresponding inhomogeneities {^} = {C,2j} U {^y-i}- The Bethe roots and 
inhomogeneities shifted by the amount of 2i are located outside. 



Here the limit in the explicit term is easy to calculate 



hm n «*-w=v™ ri(Cy-cy-i 



\<}<k<2m 



•7=1 



X EI (C2fe — C2j-l)(C2fc — C2j)(C2fc-l — C2y-l)(C2A:-l — 
l<j<k<m 

r m 

n-2isign5 7 J] fe-y 2 [fe-y 2 + 4]. (C.15) 

■;'=1 J \<j<k<m 

The combinatorial sum can be converted into a multiple integral by the same token 
as in equation (65) of [12]. We introduce a function 



N 



X (X, Q = e(C - A,)o(Q - e(A, - Q + £ - X k )w k (Q . (C.16) 



Then 

X(A.-2i,C) 



HQ 



- + 



X — t, X-^-li X — t, — A\ ^yX-Xu-^i X-Xu 



w k (Q w k (Q 



(C.17a) 



X(A. + 2i,C) 



a(C) a(C) 



+ ■ 



A,-? X-^ + 2i X-^ + 4i 



yt=l L 



X — Xk X — Xk + 4i 



(C.17b) 



Let 01 = {z £ C| 1 < |Imz| < 2}. Then 31 contains all Bethe roots. The two functions 
%(X^f2i,C,j) are meromorphic in 01. Their only poles inside 01 are all simple and are 
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located at the Bethe roots and at C,j. The corresponding residua are 

nsj^x(*-2i,Cy) = -wjfc(Cj). res x=e .x(X-2i,y = 1 , (C.18a) 

res x=Xk x& + 2iXj)=w k &j), res^.x(X + 2i,y = o(Cy) ■ (C.18b) 

Define two simple contours T*, such that (i) T + is inside the upper strip of 3? and r is 
inside the lower strip of % and (ii) all Bethe roots and all C,j, j = 1, . . .,2m, are inside 
r = r + + r (see figure 4). Decompose F in such a way that F = H* + J, where 23 
contains only Bethe roots and J contains only inhomogeneity parameters. Then we are 
very much in the same situation as in [12], and using the functions %(X ± 2i, i^) we can 
transform the right hand side of (C.14) into a multiple integral over F. As we shall see 
the notation 

gMQj^-^ { £P (C19) 
|x(© + 2i,C) ]>P 

will prove useful in that exercise. Using also (C.15) we obtain 



DP](x|Q 



n-2i<8(^-i)sign8; J] fe-y 2 [fe-y 2 + 4] 

L/=l J \<j<k<m 

r n r dco F -i N 
I II / 2;ri I I ^n(PQ) 

l i= lJ0 J V'-'V =1 ({5 + },{8-})e P2 (Z 2 ,„) 

1 2m - n card{S-}=n 



({e+},{e-})G K (Z 2m ) 



det[w/ e+ (C8+)] 

E sign(/?) g e - (eou- , £ 8 - ) . . . g - (co - , C 5 - ) jj- — - — — 



det[g (co C, ) 



2m— n 

n 

_ 7=1 1 J 



2m— n 



n ^(k+: 

7=1 ; ; 
~2m-n r dCO, 



£j <P 



£j >p 



lim 

E^0+ 



({e+},{e-})ep 2 (Z2„) 



"2m-n „ d(0 E + IT" /■ d(£> ET " 

n y 2 ^TtT ny, ni 



<j<k<2m{®j-®k-2i) 



({8+},{8-})Gp2(Z 2m ) 
card{8 }=n 



— ll/=iJr 27ti 

sign(PO) det[g e + (co E + , C 5 + )] det[g eJ (ca^ , C 8 - )] 



=det[ gj (Wj,^) 



2m— n 



2m— n 



n^ ] (co e; ) n^K + ) n^K) n^K 7 ) 

7=1 7 J L 7=1 ' J L 7=1 ' J L 7=1 1 



£~t>P 



£j<P 



E j>P 
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= lim 

e^0+ 



. U.l<j<k<2m{®j-®k-2i) 

(C.20) 



Note that we used the Laplace expansion formula for determinants in the third equation. 
To summarize up to this point, we have derived the equation 



^(x|Q % Um^n mfe _ y2 ^_ y2 



+ 41 



J-r f d(0, m, J T %Z rdft);- 

n/ r ^ [1 H[n +i / r ^ 



d<D "f(»,) 



det[g y (co y ,^]n7=i^T 



Ill<j<k<2m(®j-®k-2i) 



. (C.21) 



Here the factor YlJ=i sign 87 can be use d to reorder the columns in the determinant. 
Defining the alternating pattern 



V2*-i=£*+i(l+e), V2* = 5*-i(l+e) 



and 



we find that 



Xjk 



X(lj-2i,v k ) j = l,...,p 
%(Xj + 2i,V k ) j = p + \,...,2m, 



d e t[ gj ((O j Xk)}U si S^j = ^txY k ' e) . 



(C.22) 



(C.23) 



(C.24) 



Note that the limit £ — > 0+ is not obvious at this stage, because in the limit the 
poles of g at the inhomogeneity parameters unavoidably cross the narrow contour T. 
Below we shall widen the contour, while carefully taking account of the additional 
terms generated during this process. As we shall see, all additional terms are of order £ 
and vanish in the limit. 



C.2.2 Fusion of wave functions 

Before coming to this point we have to recall that for the spin-1 density matrix elements 
we do not exactly need (x|Q, but certain combinations of these coefficients. This 
leeds to 'fusion of the wave functions' FN, to be described in this subsection. 
Let us consider a specific matrix element 

" ! X' .' • . - ! ''p. X,„ llfn; 

(C.25) 



D[2] a 1 .... i a = Vo\T%l(M...T®£(S m )\v ) 



of the spin-1 density matrix. Since the local space is spin-1, the indices take three 
different values, otj,$k = +,0, — . The right hand side of (C.25) can be written as a 
linear combination of coefficients (x|^) which can be identified by means of (47). 

To begin with let us assume that T^i^i) = rW+^-i^W+^ + i) is contained 
in the sequence of monodromy matrix elements on the right hand side of (C.25). Then 
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we must have Xj = 2£ and Xj + \ =2t—\ for some j G {1 , . . . ,p — 1 } in all coefficients 
dI 2 1(x|^) contained in the linear combination for that specific density matrix element. 
Also 5^ = — 1 — £, and a factor 

2£— 1 2m 21-2 2m 

n(G>y-s*-2i) n nK-+i-^-2i)riK- + i-w 

£=1 /t=2f+l jfc=l fc=2^ 

= [(<oj - ^ - i) (co 7 +i - \t - i) - ie(co, - <o j+ i - 2i) + s(2 + e)] 

j+i e-i 

x EI. 11 - ^ - 3i ) - 5* - i) + e(2 + e)] 

n=j k=l 

m 

x II [(ro„-^-i)(co„-^+i) + e(2 + e)] 

k=Z+\ 

= F e X(Oj)F e ((Oj +l ) + 0(£) (C.26) 

appears. Here we used the 'spin-1 wave function' Fg defined in (65). 

In a similar way we may consider all the matrix elements of using for simplifi- 
cation the right hand side of (47). E.g. if rPiJ" (y = V2rW+(^+i)rWt(^-i) is 
contained in the sequence of monodromy matrix elements on the right hand side of 
(C.25), then 8f = 1 + £, and we have xj = 2£, Xj+\ = 21— \ for some j e{l,...,p— 1} 
and Xj = 21 for some i G {p + 1 , . . . , 2m}. Thus, there is a factor 

>/2 [(co,- - \t - i) (coy+i - ^ - i) + i(2 + e) (co,- - co,- +1 - 2i) + £(2 + e)] 
y'+i «-i 

x IT IT [ ((0 " ~ ^ ~ 3i) ((0 " - ^ - i) + e(2 + e)] 

m 

x EI [(co„-^-i)(co„-^ + i) + £(2 + £)] 

Jfc=M-l 

i-\ 

x (co; - ^ + i - ie) Yl [(co* - 5* + 3i) (co,- - ^ + i) + e(2 + e)] 

k=\ 

m 

x H [(co ! --^-i)(co / -^+i) + e(2 + e)] 

k=t+i 

= V2F e X(Oj)F f X(o j+l )F ( i(Oi) + 0(e) (C.27) 

under the integral. Again Fi and are taken from (65). We use the notation '=' 
for 'equal under the multiple integral' (C.21). The crucial point here is that the term 
proportional to (2 + e) (©/ — C0 7 +i — 2i) does not contribute under the integral (C.21) for 
symmetry considerations. This is because it multiplies a function which is symmetric 
in 0)j and co y+ i in {C.21) and the only other terms under the integral depending on co 7 
and C0 ; - +1 are det[g 7 (co 7 -, and Y\\<j<k<2m(®j ~ ®k ~ 2i). But 

K--CO;+i-2i)det[g j -(co j -,^)] 
ril<j<*:<2m(t»j-C0^-2i) 
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Matrix element Factor under the integral 





F £ ((Oy)F f ((0 7 - + i) + 0(£) 




V^F^COy^^CO^+OF^CO;) + 0(e) 




F e {Gij)F e (aj +1 )Fe{Gii)Fe(a i+1 ) + 0(e) 




V2F t (<0j) + Q(e) 




2Fe(<0j)F t (<0i) + Q(fL) 




V2F t ((0j)F e ((0i)F e ((0i + i) + 0(e) 




1 




v / 2F^(co,-) + 0(e) 




F £ (co ; -)^(co ; - +1 ) + 0(e) 



Table 1 : The polynomials under the integral. 



is antisymmetric in go, and tO/+i. 

Another example is the matrix element T^lfy) = rM+(^ - i)T^t(^e + i) + 

rW_(^ — i)rW + (^ + i) for which 8f = — 1 — e. It is the only matrix element which 
we have to express by a sum of two products of monodromy matrix elements with 
spin- j auxiliary space, and it contributes a factor 

2<?-l 2m 

n(co;-c*-2i) n 

fc=l k=2(,+l 

-2l-\ 2m 21-2 2m 

X J] (CO; - ^ + 2i) J] (CO; - Ct) + fj ~ & + 2i ) 11 ~ W 

= 2F e ((O j )F e ((O i ) + 0(e), (C.28) 

where j 6 {1, . . . ,/?} and / G {/? + 1, • ■ .,2m}. 

Proceeding case by case in a similar way we obtain the result exposed in tabular 1 . 
In the tabular it is always implied that j,j+ 1 € {1, . . . ,p} and i, i+ 1 G {p+ 1, . . . ,2m}. 
Inspecting the tabular we see that, up to corrections of the order of e, the 'wave function' 
under the integral is composed in the following way: For every zero in the sequences 
(a,-) and (0*) we obtain a factor of y/2, amounting to a total factor of 2("o( a )+MP))/2 = 
2m-n + (a)-n_(P)_ p or ever y p} us j n ( a .~j a f ac tor F( : ((Oj)Ff(oij+i) appears and for every 
zero a factor F^(o)j). From the sequence (P^) we obtain a factor F^(co,-) for every 
zero and a factor F^(co ; -)F^(co, + i) for every minus sign. This makes a total number of 
2n + (a) + n (a) +«o(P) +2«_(P) = 2m factors, and /? = 2n+(a) +n (a). With the 
factors we obtain a sequence {zj) 2 -l x , Zj £ { 1 , . . . , m} by arranging them in the order of 
ascending co y : F Zl (coi ) . . . F Zp ((O p )F Zp+l . . . F Z2m (co 2m ). 



35 



Thus, we have obtained the representation 

2 -„ + («)-n-(P) n »L 1 ^ =IJ 



[2]«i,...,a m 



ni<^< m fe-y 2 [fe-y 2 + 4 ] 



lim 

e->CH 



2m 



n 

■y=i 



dco, 



:_-wr 27ii 



det X ^' £) 



ril< 7 </t<2m((0;-(0«:- 2 i) 

[F Z1 (co x ) . . .F Zp (a p )F Zp+ , ((Op+i). . .F Z2m (co 2m ) + 0(e)] (C.29) 

for the spin-1 density matrix elements. What remains to do is to calculate the limit 
£ — > 0. For this purpose we have to deform the integration contours first. 

C.2.3 Widening the contours 

Next we would like to show that 



f2J fdaj- _ 
} = Vr27liJn 



l<y<fe<2m( ro 7 ~ (a k~ 2i) 

[F Zi (eoi) . . .^(co p )F Z; , +1 (co p+1 ) . . .F Z2m (co2 m ) + 0(e)] 

'dCO 1 r 2m 



det X J' e) 



- rdco^ - rdcoy _ 
f \ Je 2ni \ [M +l Je 2%\ \ Ui<j<k<2m(^j - - 2i) 

(C0l) . . .K ;; (co P )F Zp+1 (co p+1 ) . . .F Z2m ((o 2m ) + 0(e)] + 0(8) . 
We recall that the simple closed contours 6, S are defined in such a way that 

e = e + + e~, e = e + + e~ 



(C.30) 



e ± ,e ± cs ± 

6~ = C" , e - inside C H 
r inside 6,6. 



2i. 



2i inside 6 



(C.31) 



For C ± , C we may take large rectangles inside S 1 * 1 which are slightly narrower than 
2 in imaginary direction. The third line in (C.31) is a closed contour analogon of 
the regularization by infinitesimal shifts of the contours in [21]. The Bethe roots of 
the dominant state come in complex conjugated pairs, so-called 2-strings. For their 
enumeration we shall employ the same convention as in [21]. Those in the upper 
half plane will be labeled by odd integers and those in the lower half plane by even 
integers. By definition the contours 6 and 6 encircle not only all Bethe roots Xj and 
all inhomogeneities Vj but also the down-shifted upper Bethe roots X,2 7 -i — 2i and the 
up-shifted lower Bethe roots + 2i as well as the down-shifted upper inhomogeneities 
V2y-i — 2i and the up-shifted lower inhomogeneities V2j + 2i. 
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In preparation of the proof of (C.30) we introduce the notation 



/(©i,..., (02m) = 



F ZI (coi) . . .F Zp (n p )F Zp+l (co p+1 ) . ■ .F Z2m (co 2m ) + 0(e) 

Ill<j<ifc<2m( (B y- (B Jt-2i) 



(C.32) 



where the polynomial, including the (e) contribution, is the same as under the integrals 
in (C.30). Then the integral on the right hand side of (C.30) can be written as 



We shall show that, if we successively replace the integrals in this expression by integrals 
over r, the total error will be of the order £. 

(a) For the rightmost integral we note that / considered as a function of iO p+ i is 
holomorphic inside 6. There is a factor (C0i — (O p +i — 2i) . . . (<O p — co p+ i — 2i)((O p+ \ — 
— 2i) . . . (cOp+i — (02 m — 2i) in the denominator, but with our choice of contours 
(Oj — 2i is outside C for (Oj EG, j = l,...,p, and the same is true for (Oj + 2i for (Oj E C, 
j = p + 2, . . .,2m. The function %(iO p+ i + 2i,Vq( p+ u) has outside F but inside C at 
most a single pole occurring at Vq( p +i) — 2i if Q(p + 1) is odd. Then there is I E Z m 
such that Vj2( p _|_i) — 2i = ^ — i + is and, if we contract the contour from C to T, the 
pole contributes a term having a factor F Zp+l (^ — i + is) = 0(e) in the numerator. 
Hence, the numerator of this term is 0(e). In the denominator we have factors of 
ro i ~ v G(p+i) = ro 7 -i-iefor j = 1, . . . ,p or v Q(p+1) - CO; -4i = ^ - (Oj - 3i + is 
for j = p + 2, . . . , 2m. It follows that the absolute value of the denominator is bounded 
from below if the co ; are on C for j = 1 , . . . , p or on C for j = p + 2, . . . , 2m. Hence, 
the additional term that may be generated by contracting the contour from C to F will 
at most contribute to order s, even after performing the summation and the remaining 
integrations in (C.33). We symbolize this by writing 




(C.33) 




X(co p+ i +2i,v G(p+1) )/(coi, . . . ,co 2m ) . (C.34) 



— E 



(b) In order to proceed by induction we define 



k =/(coi, ...,co 2m ), 



(C.35a) 




(C.35b) 



We want to show that 




(C.36) 
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We have already shown that this is valid for n = 1 . To proceed further we have to know 
the structure of /„. 

(c) As a preparatory step let us consider the left hand side (C.36) for n = 2. Then 



Here we have used that / is holomorphic as a function of (Q p+ i for (O p+ \ on and 
inside F. We insert 1\ into the left hand side of (C.36) for n = 2 and contract the 
integration contour from S to T. Due to our special choice of the contours S, S 
the poles at C0 p + 2 = COy — 2i, j = 1, . . . ,p and at iO p+ 2 = CO/ + 2i, j = p + 3, . . .,2m, 
remain outside the contour during the process of deformation. The only pole of 
X(g>p+2 + 2i, Vq( p+ 2)) which may be crossed is, in case that Q(p + 2) is odd, a simple 
pole at Vg( p+2 ) — 2i. The situation is the same as in case (a) above. And as above we 
can see that such a term gives only an order-£ contribution, even after performing the 
sum and the remaining integrals: Vg( p+2 ) = ^ + i + ie => F z 2 (Vg( p+2 ) — 2i) = (e) 
and in the denominator we may have (Oj — Vg( p+2 ) = 00/ — \i — i — is for j = 1 , . . . , p or 
Vq(p+2) ~ G>7 - 4i = ^ - C0 7 - - 3i +is for j = p + 3 , . . . , 2m, as before, or X kp+l - V G(p+2 ) 
or v 2(p+i) ~~ v e(p+2)- The latter terms are of no danger, since we assume that the £,j are 
mutually distinct and distinct from all Bethe roots. Thus, we see that the same argument 
as above works. 

However, we now have singularities of 1\ which are crossed in the course of the 
deformation of the contour. They give additional contributions. The summand with 
k p+ \ = j has a term Xj — (£) p+ 2 — 2i in the denominator, giving rise to a simple pole in 
iO p+ 2 at Xj — 2i, if j is odd. When contracting the contour of the a)p + 2-integral this pole 
causes a contribution (see (C.16), (C.8)) proportional to 



which is symmetric in p+ 1, p + 2. Such a term vanishes under the sum over all 
permutations in (C.33). Another additional contribution comes from the second term in 
(C.37), which has a factor Vgi p+ u — iO p+ 2 — 2i in the denominator. So we have a pole 
at Vg(p + i) — 2i. It will give only 0(e) -corrections when we integrate over (O p+ 2, as we 
have already seen. 

(d) In the general case the argument is very similar. 7„_i is obtained by iterating the 
integrations over F. In every integration the 7 7 _ i under the integral is holomorphic on 
and inside T by construction. Hence, we obtain a sum over the pole contributions of 
X((O p+ j + 2i,VQ(p + y)). This means that 7 ra _i is a sum over terms of the form 

t(®p+n) = /(COi ,...,(O p ,X p+ i,... ,X p+n -i, C0p+„, . . . , (02m) y P +l ■ ■ ■ 1 , (C.39) 

where xj is either a Bethe root or Vqj. If Xj is a Bethe root, say kg, then the corresponding 
factor yj = we(\Qj), if xj = Vqj, then yj = a(vgy). Let us consider 



N 



k p +i=l 



+/(coi , . . . , (a p , Vq(j,+i) , (Op+2, • • • , co 2m )a(v G (p +1 )) . (C.37) 



X&j,VQ(p+2))wj{V Q{p+l) ) = d{Xj)wj{v Q{p+ 2))wj(VQ (p+l) ) 



(C.38) 




(C.40) 
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If we shrink the contour to T, we obtain at most one pole contribution from % which is 
at Vg( p+ „) — 2i if Q(p + n) is odd. This contribution is 0(e) by the same argument as 
above. Also for the contributions stemming from t we can argue as above. Note that we 
do not have to consider double poles (or poles of even higher order), since, if Xj and 
Xk are the same Bethe root, Tig say, then t((0 p+ „) has a factor wg(VQj)wg(VQk), which 
vanishes under the antisymmetrizing sum in (C.33). Thus, we have accomplished the 
proof of (C.36). 

(e) It follows from (C.36) that we can replace all the C-integrals in (C.33) by T- 
integrals. The total error will be only of order e. To finish the proof of (C.30) we have 
to proceed with the contours C. Let us set 

J p+ i=hm- P , (C.41a) 
r dco„ . 

Jn = %(co„-2i,v Q „)7„ + i, n = /?,..., 1. (C.41b) 
We want to show that 

f ^X(»« " 2i,v e „)/ n+ i = e J n . (C.42) 

The proof is very similar as before. 

(f) Let us start with n = p. Then J p+ \ = hm-p is a sum over terms of the form 

s (®p) =f(®l,---,®p,Xp+l,---,X2m)y P +l---y2m, (C.43) 

as follows from (C.39). If we shrink the contour in 

f e ^X{v>p-X,VQpHa>p) (C44) 

to r, we obtain at most a single pole contribution from % stemming from a pole at 
Vq p + 2i in case that Qp is even. Then Vq p + 2i = ^ + i — ie for some £ € Z m . But 
now F z + i — ie) = 0(e), and the numerator in the generated term is 0(e). The 
denominator contains terms co ; — (vq p + 2i) — 2i = (Oj — \g — 3i + ie which are bounded 
from below in the absolute value for CO/ G C, j = 1, . . . ,p— 1. Hence, the whole term is 
of order e, even after integration and summation, and can safely be forgotten. As for 
the singularities of 5(co p ) we have again two types. If xj is a Bethe root, say Tig, then a 
term (O p — Xe — 2i occurs in the denominator. It comes together with a factor wg(VQj). 
When calculating the residue at Tig + 2i, which is non-zero only if i is even, we obtain 
something proportional to 

X(X e ,VQ P )wi(v Qj ) = a'(Xi)w e (vQ P )w e (v Qj ) (C.45) 

which yields a term that vanishes under the sum in (C.33) due to symmetry reasons. 
Double poles can be excluded by the same argument as above. Finally we may have 
Xj = Vqj. Then a factor co^, — Vqj — 2i is present in the denominator resulting again at 
most in an 0(e) -contribution. 
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(g) Iterating the above arguments we obtain (C.42), and the proof of (C.30) is 
complete. With (C.30) we can now perform the limit £ — > 0+, because it is trivial at the 
right hand side of the equation. With the definition 



we obtain 



lim 



2m 

n 



d(0; 



yiP) 



det X J £) 



E->0+ 



(C.46) 



ni<y</K2m(<»/-<0*:- 2 i) 

[F Z1 (a)!) . . .F Zp ((o p )F Zp+l (a> p+1 ) . . .F Z2m (co 2m ) + 0(s)] 
do,- , J r ^ /■ dco; 



" p /■ dco,- , J T 3^ r dco,— , 



detxj 



J Ul<j<k<2m(^j - COfc - 2i) ' 

With this we have derived the multiple integral representation 

2~n+ («)-"- (P) 



(C.47) 



[2]ai,...,a m 



ni<;<.< m fe-y 2 [fe-y 2 +4] 



t^t /" d(0/ . . f dCO;— , 

FT / irrFzA^i) 1 I / ^ F z, (»/ 



: '" '• do. 



(p) 



detx}; ; IT/ 



=1 



for the inhomogeneous density matrix of the isotropic spin-1 chain. 



(C.48) 



C.3 The linear integral equations 

In this subsection we shall derive a pair of coupled integral equations for the functions 
%(k± 2i, v). For this purpose we first of all note that the function % has been defined in 
(C.16) in such a way that 

X(^,v) = c'(A,> y (v), j=l,...,N (C.49) 

(see (C.8)). We shall use some of the functions that appear in the formulation of the 
thermodynamics of the model [30] and that are collected in appendix A. In first place 
we need the functions 21 and 21 (see (A.7a)). It follows from their definition and from 
the fact that AM has no zeros in S + U S~ that the only zeros of 2l(X, + 2i) in S~ are 
simple zeros at X,2 7 _i — 2i, while the only poles in S~ are simple poles at 7i2j. Similarly, 
the only zeros of 2l(X — 2i) in § + are simple zeros at %2j + 2i and its only poles in § + 
are simple and located at Xij- 1 . 

Hence, using (C.17b), (C.49), we conclude that the only singularities of the function 
X(X + 2i,^+)/2l(X + 2i) inside S" are 

(i) a simple pole at X = with residue — 1 , 

(ii) simple poles at X = 7i2k-i — 2i with residua W2k-i (^ + )- 
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Similarly, all singularities of %(X — 2i,^ + )/2l(X — 2i) inside § + are (see (C.17a) and 
(C.49)) 

(i) a simple pole at X = ^ + with residue 1 /2t(^~), 

(ii) simple poles at X = Xn + 2i with residua — W2/t(^ + ). 
With this it follows by means of (C.17) that 

1 aOr) 1 



x-^+-4i 

r ^X^2i^) r d^ ^+2ij^) 

7e+ 2jri j%-2i) V ; h~ 2%\ a(/i + 2i) v M ; ' 

S(^) 1 oft") 1 



2l(^")X-^- 2l(^")X-^-+4i 

_,d^-2^+) K{X )+ f ^5fc2i^) 
7e + 27ii a(/i-2i) V 7 7e-27ti 2l(/i + 2i) v w 

for X G C + in the first equation and X G C~ in the second equation. 

If now v = ^ we have to repeat almost the same considerations. The function 
X(X + 2i,^ _ )/2l(X + 2i) has the following singularities inside S~: 

(i) an simple pole at X = with residue ct(^~)/2l(^ + ), 

(ii) simple poles at X = A^-i — 2i with residues W2k-i (^~)- 
Furthermore, all singularities of %(X — 2i, J;~)/2l(A, — 2i) inside S + are 

(i) a simple pole at ^ + with residue — a(£~), 

(ii) simple poles at X = X^k + 2i with residua — W2k(£,~)- 
Thus, we obtain 



2l(^+)X-^+-4i K(Z+)X-Z+ x-t,- 
7e+ 27ti j%-2i) V HJ Je-2ni 21(^ + 21) v H h 

a{%-) i i , sir) 



21(5+) A, -5+ A,-£-+2i X-^-+A\ 

_ [ ± X(M-2i,%-) T ^X|+2i^) 

7e + 27ii a(/i-2i) ye- 2m 2l(/i + 2i) v 

where X G C + in the first equation and X G C~ in the second equation. 
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After taking the Trotter limit N — > °° the driving terms in the integral equations 
above are singular at ^ = 0. The singularity can be removed by taking appropriate linear 
combinations. We define 

g±(uA = (m-)im-) aft + )/»ft-A to=F2i,§ + A 

Then, using also (A.9a)-(A.9c), we arrive at equations (66), (67) of the main body of 
the text. 

In order to express the determinant under the integral in (C.48) in terms of the 
functions G ± , 5 ,± we use the matrix &^ defined in (69). Because the determinant of 
the matrix in (C.52) is 2/58(£ - ), we obtain 

detxjj? • (C53) 
Inserting this into (C.48) we arrive at the multiple integral representation (70). 



det© 



(p) 



2i 



i n i 
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